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ABSTRACT 

A  numerical  approximation  scheme  for  the  estimation  of  functional  parameters  in 
Euler-Bemoulli  models  for  the  tansverse  vibration  of  flexible  beams  with  tip  bodies  is  developed. 
The  method  permits  the  identification  of  spatially  varying  flexural  stiffness  and  Voigt-Kelvin 
viscoelastic  damping  coefficients  which  appear  in  the  hybrid  system  of  ordinary  and  partial 
differential  equations  and  boundary  conditions  describing  the  dynamics  of  such  structures.  An 
inverse  problem  is  formulated  as  a  least  squares  fit  to  data  subject  to  constraints  in  the  form  of  a 
vector  system  of  abstract  first  order  evolution  equations.  Spline-based  finite  element 
approximations  are  used  to  finite  dimensional ize  the  problem.  Theoretical  convergence  results  are 
given  and  numerical  studies  carried  out  on  both  conventional  (serial)  and  vector  computers  are 
discussed. 


1.  Introduction 

We  develop  here  numerical  approximation  methods  for  the  estimation  of  functional  or  more 
precisely,  spatially  varying  parameters  that  describe  material  properties  in  continuum  models  for 
elastic  structures.  In  particular,  we  consider  the  identification  of  the  flexural  stiffness  and 
Voigt-Kelvin  viscoelastic  damping  coefficients  in  Euler-Bemoulli  models  for  the  transverse 
vibration  of  long,  slender,  flexible  beams  with  tip  appendages.  The  primary  motivation  for  the 
study  we  report  on  here  is  the  modeling  and  ultimately  the  control  of  the  dynamics  of  large  flexible 
spacecraft  The  type  of  structures  to  which  we  are  referring  includes  satellites  with  flexible 
appendages  (solar  panels  and  the  like)  antennas  (reflectors  as  well  as  supporting  structures)  and 
trussed  masts  and  platforms,  both  shuttle  attached  and  free  flying. 

The  difficulties  involved  in  die  design  of  efficient  and  practical  control  laws  and  in  particular  the 
need  for  extremely  high  fidelity  models  for  structures  of  these  types  are  well  documented  (see,  for 
example,  [1],  [8],  [21],  [22]).  Their  high  flexibility,  light  damping,  construction  with  new  and 
relatively  untested  composite  materials  (usually  graphite-epoxy)  and  overall  complexity  together 
with  their  use  in  a  fuel  limited  and  highly  variable  environment  all  contribute  to  making  space 
structure  stabilization  and  control  a  formidable  task.  It  is  becoming  increasingly  clear  that  the  use 
of  continuum  or  distributed  models  with  spatially  and  /  or  temporally  varying  functional  parameters 
has  the  potential  to  offer  several  distinct  and  significant  advantages.  Included  among  them  is  the 
ability  to,  in  some  sense,  capture  the  physics  and  inherent  infinite  dimensionality  of  the  dynamics 
while  at  the  same  time  greatly  reducing  the  number  of  unknown  or  experimentally  indeterminable 
material  parameters  which  have  to  be  identified  (see  [15],  [18],  [23],  [28],  [35]). 

In  our  study  we  have  considered  exclusively  Voigt-Kelvin  viscoelastic  damping  which  is  based 
on  the  hypothesis  that  the  damping  moment  is  proportional  to  strain  rate.  There  exists  considerable 
evidence  to  suggest  that  damping  mechanisms  in  composite  materials  are  significantly  more 
complex  than  the  one  described  by  the  Voigt-Kelvin  model.  For  example,  it  has  been  conjectured 
by  some  investigators  that  an  appropriate  model  might  involve  hysteretic  or  hereditary  effects. 
However,  since  there  are  a  number  of  materials  for  which  the  Voigt-Kelvin  assumption  is 
appropriate  and  moreover,  since  at  present  many  questions  regarding  the  modeling  of  structural 


damping  mechanisms  remain  open,  we  feel  that  the  Voigt-Kelvin  model  leads  to  a  reasonable  class 
of  examples  and  problems  on  which  we  can  begin  to  develop,  test,  and  evaluate  identification 
schemes. 

Our  treatment  here  is  similar  in  spirit  to  some  of  our  earlier  efforts  and  the  work  of  others  on 
inverse  problems  for  elastic  structures  (see  [2],  [3],  [4],[5],  [6],  [14],  [17],  [26],  [31]). 
Formulating  the  identification  problem  as  a  least  squares  fit  to  data,  the  scheme  we  develop 
involves  a  spline  based  finite  element  approximation  to  the  hybrid  system  of  coupled  ordinary  and 
partial  differential  equations  describing  the  dynamics  of  the  structure  together  with  a  spline  based 
discretization  of  the  admissible  parameter  set 

Our  approach  here  specifically  differs  from  the  one  taken  in  [5],  [6]  in  that  the  present  scheme  is 
derived  from  an  alternative  state  space  formulation  for  the  underlying  dynamical  equations.  We 
consider  the  higher  order  analog  of  the  classical  conservative  formulation  for  a  second  order 
hyperbolic  equation  as  a  first  order  vector  system  in  the  natural  states  of  strain  ux  and  velocity  ut. 
We  have  considered  identification  schemes  based  upon  this  formulation  previously  in  [31]. 
However  by  replacing  the  semigroup  theoretic  convergence  arguments  used  there  with  weak  or 
variational  arguments  (in  the  spirit  of  those  commonly  found  in  the  finite  element  literature)  as  used 
in  [5],  we  are  able  to  significandy  weaken  the  hypotheses  necessary  to  ensure  convergence.  We 
point  out  below  that  the  weakening  of  these  hypotheses  has  both  theoretical  and  computational 
significance. 

Along  with  reporting  theoretical  convergence  results,  we  discuss  numerical  findings.  Our 
computational  results  are  based  upon  extensive  numerical  studies  which  involved  a  variety  of 
examples  and  two  machines.  In  addition  to  testing  our  scheme  on  a  conventional  serial  computer 
(an  IBM  3081)  we  vectorized  our  codes  for  the  Cray  1-S  and  then  benchmarked  some  of  our  runs 
in  order  to  explore  the  potential  of  vector  architectures  in  the  context  of  inverse  problems  for 
systems  described  by  distributed  parameter  models. 

We  provide  a  brief  outline  and  summary  of  the  remainder  of  the  paper.  In  Section  2  we  specify 
the  ordinary  and  partial  differential  equations  which  govern  the  underlying  dynamics  of  the 
structure  and  precisely  formulate  the  identification  problem.  We  reformulate  the  initial-boundary 
value  problem  as  an  abstract  second  order  evolution  equation  and  then  as  a  first  order  vector 
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system.  Existence,  uniqueness  and  regularity  results  for  solutions  are  summarized.  Section  3 
contains  the  abstract  approximation  theory  and  convergence  results.  A  spline-based  scheme  is 
discussed  in  detail  in  Section  4  and  our  numerical  findings  are  reported  and  summarized  in  Section 
5. 

We  use  standard  notation  throughout.  For  X  and  Y  Banach  spaces,  the  Banach  space  of 
continuous  linear  transformations  from  X  into  Y  is  denoted  by  £  (X,Y).  When  X  =  Y  we  use  the 
shorthand  notation  £  (X).  The  spaces  of  (equivalence  classes  of)  functions  f  from  an  interval  4 
into  X  which  satisfy 

I  f(6)  |  x  d0  <  oo  or  ess  sup  If(0)lx<~ 

d 

are  denoted  respectively  by  L^d  ;  X)  and  L^d  ;  X).  For  k  =  0,1,2...  the  space  of  X-valued 
functions  with  k  continuous  derivatives  on  d  are  denoted  by  Ck(d  ;  X).  When  k  =  0  we  use 
C(d  ;  X).  The  completion  of  the  space  C^d  ;  X)  with  respect  to  the  norm 

\  j_0  t 

is  denoted  by  Hk(d;  X).  When  X  =  R  we  use  simply 

2.  The  Identification  Problem 

We  consider  the  identification,  or  estimation,  of  the  mass  and/or  material  properties  of  a  long, 
slender,  flexible,  viscoelastic  beam  of  length  l  and  spatially  varying  mass  density  p  which  is 
clamped  at  one  end  and  free  at  the  other  with  a  body  rigidly  attached  at  the  free  end  (see  Figure  2. 1 
below). 


1^(3),  Uo(d),Ck(d)  and  Hk(d). 
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Figure  2.1 


We  assume  that  the  material  behavior  of  the  beam  is  that  of  an  idealized  Voigt-Kelvin  solid  with 
modulus  of  elasticity  E  and  coefficient  of  viscosity  CD  (see  [30  ]).  We  assume  further  that  E,  CD 
and  the  cross  sectional  moment  of  inertia  I  of  the  beam  are  in  general  spatially  varying.  We  take  the 
mass  properties  of  the  tip  body  to  be  mass  m  and  moment  of  inertia  J  about  the  center  of  mass  O 
which  is  assumed  to  be  located  at  a  distance  c  from  the  tip  of  the  beam  directed  along  the  beam's  tip 
tangent  (see  Figure  2.2  below). 


=  .,©‘ 


Figure  2.2 


We  note  that  there  is  no  essential  loss  of  generality  in  assuming  that  the  mass  center  of  the  tip  body 
is  not  offset  from  the  dp  tangent  of  the  beam.  We  refer  the  interested  reader  to  [31]  where  the 
more  general  situation  is  treated.  Also,  the  problem  with  non-zero  mass  center  offset  can  be 
transformed  into  a  problem  of  the  general  form  of  the  one  which  will  be  considered  here.  See  [32] 
for  details. 

Letting  u  =  u(t,x)  denote  the  transverse  displacement  of  the  beam  at  position  x  at  time  t  and 

assuming  only  small  deformations  ( |  u(t,x)  |  «  l,  |  —  (tpc)  |  «  1),  the  Euler-Bemoulli  theory 

dx 

and  elementary  Newtonian  mechanics  yield  the  hybrid  system  of  ordinary  and  partial  differential 
equations  (see  [19],  [34]) 


a2  2  -.2  3 

d  u  d  3  u  d  u 

(2. 1 )  p— xr(t,x)  +  {El— (t,x)  +  CDI—— -(t,x) }  = 

dt  dx  3x  dx  dt 


3  du 

T~  O— <t,x)  +  f(t,x), 
dx  dx 


xe(0,-t),  t  >  0 


(2.2)  m^(t,4  )  +  me  (t ,Jt)  -  -2-  (Elli  + 

dt  dt  dx  dx 


d3u 


,d2u 


d3u  du 

CDI— -)(t,i)  =  o— (t,i)  +  g(t), 

dx2dt  3x 


t>0 


(2.3)  mci^(t,i)  +  (J  +  me2)  (t  ,4)  +  EI^  M)  + 

dt  dt  dx  dx 


d  u  du 

Cy-V-(t,4)  =  -CO— (t,i)  +  h(t), 

dx  dt  dx 


t>0 


(2.4)  u(t,0)  =  0,  ^<t,0)  =  0,  t  >  0 

dx 
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(2.5)  u(0,x)  =  <|>(x). 


—  (0,x)  =  v(x), 
dt 


x  e  [0,£]. 


Equation  (2.2)  and  (2.3)  are  derived  from  the  usual  transverse  and  rotational  equilibrium 
considerations  at  the  free  end.  The  geometric  boundary  conditions  (2.4)  are  the  zero  displacement 
and  zero  slope  constraints  at  the  clamped  end.  The  functions  f  =  f(t,x),  g  =  g(t),  h  =  h(t)  and 
a  =  cr(t,x)  denote  externally  applied  loads  in  the  form  of  moments  (h)  and  transversally  (f  and  g)  or 
axially  (a)  directed  forces  exerted  on  the  beam  or  tip  body.  (In  fact,  h(t)  =  h(t)  +  cg(t)  where  h  is 
an  externally  applied  torque  on  the  tip  body).  The  temporal  boundary  conditions  (2.5)  reflect  the 
initial  displacement  and  velocity  distributions  which  are  assumed  to  be  given  by  the  functions  and 
y  respectively. 

We  treat  the  initial-boundary  value  problem  (2.1)-  (2.5)  in  the  form  of  an  abstract  second  order 
evolution  equation  which  we  then  rewrite  as  an  equivalent  first  order  vector  system.  The  particular 
state  space  formulation  we  choose  forms  the  basis  for  the  finite  dimensional  approximation 
schemes  we  develop  in  the  next  section.  It  also  allows  us  to  easily  establish  existence,  uniqueness 
and  regularity  results  for  solutions  to  (2.1)  -  (2.5)  using  the  theory  of  abstract  parabolic  systems. 

Let  H  denote  the  Hilbert  space  R2  x  1^(0, -6)  with  inner  product 

(Tl2^2*e2>H  =  Wl  +  5i^2  +  <M2>0 
and  let  V  denote  the  Hilbert  space 

V  =  {(ti,5,0)  e  H :  6  e  H2(0,*),  0(0)  =  D0(O)  =  0,  -q  =  0(^),  $  =  D0(*)} 
with  inner  product 

<^j,  ^2>v  =  <EI(D20i),  D202>o 

for  Oj  =  (0j(£),  D0j  (-6),  0j)  e  V,  i  =  1,2.  In  the  above  definitions  the  inner  product  <-,->0  is 
the  standard  one  on  L2(0,  l)  andD  denotes  the  spatial  differentiation  operators  or  With  H  as 

the  pivot  space,  we  obtain  the  usual  dense  embeddings  VcH  =  H'cV'. 

We  consider  the  system  (2.1)  -  (2.5)  in  the  form  of  the  abstract  second  order  initial  value 


problem 

(2.6)  7H0  utt(t)  +  C0ut(t)  +  Xqu(i)  =  B0(t)u(t)  +  O^t),  t  >  0 

(2.7)  u(0)  =  $,  ut(0)  =  v 

in  the  state  u(t)  =  (u(t,-fc),Du(t,l),u(t,  •))  e  H.  The  abstract  mass,  damping  and  stiffness 

operators  TUq,  and  OCq  are  given  formally  by 

TTUo(ri,^,0)  =  (mrj  +  me!;,  merj  +  (J  +  mc^)!;,  p0) 

G0(tU,6)  =  (-D(CdI(D20))(  t ),  CdI(D20)(^),D2(CdI(D20))) 
and 

%o(n,te)  =  (~d(ei(d20))(  l ),  ei(d20)(2),  d2(ei(d20))) 

respectively.  For  each  t  >  0,  the  operator  valued  function  B0  and  input  or  forcing  function  STq 
take  on  the  values 

Bo(t)(Ti£,0)  =  (-o(t,-t)(D0(*)),  -  c<J(t,-e)(D0(*)),  D(a(t,)(D0))) 
and 

3ro(t)  =  (g(t),h(t),f(t,0). 

The  initial  conditions  $  and  \j/  are  given  by 

$  =  «K*),D<K*),<|>) 

and 

The  formal  definitions  given  above  can  be  made  precise  and  the  existence  and  uniqueness  of 
solutions  to  the  initial  value  problem  (2.6),  (2.7)  can  be  established  if  we  make  the  following 
assumptions. 

Aj  The  functions  p,  El  and  CDI  are  elements  in  C[0,-t]  and  there  exists  a  positive  constant  a  for 
which  p(x)  £  a,  EI(x)  £  a,  CDI(x)  £  a,  x  e  [0,  ■£]. 
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Aj  The  mapping  t  -*  o(t,)  is  an  element  in  ^((O.T);  H^O.-t))  for  some  T  >  0. 

A3  The  mapping  t  ->  f(t,)  is  an  element  in  L^ftO.T);  1^(0, Z))  and  g,h  e  L^O.T). 

A4  The  function  <j>  is  an  element  in  H2(0,  Z )  with  <|>(0)  =  D0(O)  =  0  and  y  e  1^(0,  l )  with  y(  l ) 
and  Dy(-t)  defined. 

Under  the  hypotheses  A1  -  A4  above,  the  operator  TUg  is  a  bounded  linear  operator  from  H  onto  H 
andG0:  Dom  (Gq)  c  H  — >  H  and  OCq:  Dom  (OCq)  c  H  — *  H  are  densely  defined,  nonnegative, 
self  -adjoint  operators  defined  on  Dom(Co)  =  {6  e  V  :  CDI(D20)  e  H2(0,t)}  and  Dom  (OCq)  =  { Q 
e  V:  EI(D20)  e  H2(0,l)}  respectively  (see  [32]).  For  each  t  e  (0,T),  fa0(t)  e  £(V,H)  and  D^t)  e 
H  while  $  e  V  and  y  e  H.  It  also  follows  that  fag  e  L^ftO.T);  23(V,H))  and  STg  e  L2((0,T);  H). 
We  shall  call  a  mapping  t  ->  u(t)  from  [0,T]  into  H  a  strong  solution  to  (2.6),  (2.7)  if 

u  e  C  ([0,T];  V)  n  C^O.T};  V)  n  Cl([Q,T]\  H)  n  C2((0,T];  H), 

u(t)  e  Dom  OCq),  ut(t)  e  Dom  (T5 q),  t  e  (0,T|,  and  u  satisfies  (2.6)  and  (2.7)  where  the  time 
derivatives  are  interpreted  in  a  strong  (norm)  sense  in  H.  We  shall  call  a  mapping  t  -»  u(t)  from 
[0,T]  into  H  a  weak  solution  to  (2.6),  (2.7)  if 

u  e  C([0,T];  V)  n  H^O.T);  V)  n  C^tO.T];  H)  n  H2((0,T);  V') 

and  it  satisfies  the  initial  value  problem  (2.6),  (2.7)  with  the  operators  Cg  and  DCg  replaced  by  their 

natural  extensions  to  operators  in  X(V,V')  and  the  time  derivatives  are  interpreted  in  a  weak  or 
distributional  sense  (see  [20  ],  [27]).  A  function  u  =  u(t,x)  will  be  called  a  strong  (weak)  solution 
to  the  initial-boundary  value  problem  (2.1)  -  (2.5)  if  the  mapping  t  -»  u(t)  given  by 
u(t)  =  (u(t,£),Du(t,4),  u(t,  ))  is  a  strong  (weak)  solution  to  (2.6),  (2.7). 

Our  approximation  theory  for  the  estimation  problem  to  be  developed  below  is  based  upon  the 

reformulation  of  the  initial  value  problem  (2.6),  (2.7)  as  a  first  order  vector  system.  This 
reformulation  is  formally  equivalent  to  rewriting  the  initial-boundary  value  problem  (2.1)  -  (2.5) 


as  a  first  order  system  in  the  states  iAi  (strain)  and  ut  (velocity)  (see  [3],  [31]).  We  note  that 

since  the  stiffness  operator  OCq  is  nonnegative  and  self-adjoint  it  has  a  unique  nonnegadve, 

1/2 

selfadjoint  square  root  OCq  :  V  c  H  — >  H.  It  can  be  written  in  factored  form  as 

OCq  =  L|jL 

where  L:VcH->  1^(0,  4)  is  given  by 

L$  =  Dty 

for  6  =  (0(4),D0(4),  6)  e  V,  and  L|j :  Dom  (L|,)  c  ^(0,4)  -»  Hby 
Dom(l4)=  (ee^O, 4)  :  EI0  e  H^O.t )} 


(2.8)  L|j0  =  (-D(EI0)(4),  ElQil),  d2(EI0)). 

If,  for  t  e  C[0,4]  with  x(x)  £  a  >  0,  x  e  [0,4],  we  let  L2tX  denote  the  Hilbert  space  L2(0,4) 
endowed  with  the  inner  product 

<®I»®2>0,T  “  <‘c®i»®2>0 

then  L*  given  by  (2.8)  with  El  replaced  by  x  is  the  Hilbert  space  adjoint  of  L  as  a  mapping  from 
V  c  H  into  . 

We  note  that  L  e  XfV.L^gj)  is  a  Hilbert  space  isomorphism  with 


<^i’  ^2>v  *  <*o  *o  =  <L§1,L§2>0.ei 


and  L'* :  1^(0, 4)  — »  V  given  by 


*  *  2  .  * 

L_10  =  (J  J  0(y)dydx,  j0(x)dx,  J  J  0(y)dydx). 

oo  o  oo 


We  also  have 


Letting  %  ■  Lj(0 ,t)  x  H  with  inner  product 

(2.9)  <(01,(lli  4l.Xl)).(02»(1l2^»X2))>H  *  +  (n2»42»X2)>H 

and  V  =  LjCO.i)  x  V  with  inner  product 

<(® i »X i)»(®2* X2^ V  *  El +  <Xi»X2>v 

we  have  the  dense  imbeddings  VcKcV.  We  consider  the  initial  value  problem  for 

z(t)  =  (w(t),  v(t))  £%  given  by 

(2.10)  wt(t)  =  Lv(t) 

(2.11)  m0  v(t)  =  - 14  w(t)  -  L* ^(t)  +  (t)L'1  w(t)  +  S' Q(t)  0 < t  £ T 

(2.12)  w(0)  =  L$,  v(0)  =  V 

which  we  rewrite  as 

(2. 13)  ZjO)  =  Cf(t)z(t)  +  y (t),  0  <  t  £  T, 

(2.14)  z(0)  =  Zq 
where  . 

(2.15)  C|(t)  -  Vi  + 13  (t) 

with  Ql :  «©  c  %  -rtC,  lie  L^O.T);  33(H)),  y  e  L2((0,T);  %)  and  Zq  e  %  given  by 

Cl(©»X)  =  (Lx>  “  Lgj0  ~  ^»o 
for  (0,x)  e  JS  »  Dom(Lg!)x  DomfCg), 

W(t)(0,(ri,4,X))  =  (0,Hlo  B^OL'1©), 

for  (0,(Tl^,x))  £  %, 

y(o  =  (o,mo  y0(o) 

and 

Zq  =  (L$,y) . 

In  formulating  the  inverse  problem,  we  keep  technical  details  to  a  minimum  by  considering  only 
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the  estimation  of  the  beam's  spatially  varying  flexural  stiffness  El  and  viscous  damping  coefficient 
CDI.  Extending  the  finite  dimensional  approximation  methods  and  corresponding  convergence 
theory  which  are  developed  below  so  as  to  be  applicable  to  the  identification  of  other  structural  or 
input  parameters,  for  example  mass  properties  (of  the  beam  and/or  tip  body),  initial  conditions  or 
loading,  is,  at  least  in  principle,  routine  (see  [7]  [  12]  [14]  [16  ]  [  31]). 

Let  =  C[0,t]  x  C[0,l]  with  norm 

(2.i6)  kl<\  =  l(qi^2)l<vas  lqJ~+  lq2L 

=  sup  lq.(x)|  +  sup  | q2(x) | . 
ttfO.i]  ufo.l) 

We  take  the  admissible  parameter  space  Q  to  be  a  compact  subset  of  (compact  with  respect  to  the 
metric  topology  induced  by  the  norm  (2.16)).  Recalling  assumption  (i)  we  assume  further  that  the 
set  Q  has  the  property  that  all  q  =  (qt,q2)  £  Q  satisfy  qj(x)£a  and  q2(x)  £  a ,  x  e  [0,4]. 

We  formulate  the  identification  problem  as  a  least-squares  fit-to-data  over  the  admissible  parameter 
space  Q.  We  assume  that  the  structure  has  undergone  a  time  varying  elastic  deformation  in 
response  to  the  initial  conditions  described  by  $  and  y  and  the  input  loads  represented  by  f,g,h  and 

o.  Denoting  the  observation  space  by  Z,  we  assume  that  at  times  tj,  i  *  1,2 . v  measurements 

£  (tj)  e  Z  (e.g.  displacement,  velocity,  slope,  strain,  etc.)  were  taken  from  the  structure. 

We  require  that  Z  be  a  linear  space  endowed  with  a  norm  |  |  ^  and  let  T  denote  an 
appropriately  defined  continuous  mapping  from  K  into  Z  .  For  example,  suppose  that 
displacement  measurements  have  been  taken  at  the  points  xj,  j  *  1,2,.. .41  along  the  span  of  the 
beam.  We  choose  Z  as  Euclidean  p-space,  *  and  take  T  to  be 

r(z)  =  (e(x,),8(x2) . 8(x^))T 

where  z  =  (w,v)  £  %  and 

©  =  (8(4),  D8(t),8)  =  L*1weV. 

With  distributed  strain  or  velocity  observations,  we  would  take  T(z)  =  w  or  T(z)  =  v  respectively. 
We  formulate  the  identification  problem  as  follows 


(ID)  Given  C(tj)  e  Z,  i  *  find  q  e  Q  which  minimizes 

*(q)  =  £  lr(z(t.;q))-C(t)|2z 

where  for  each  q  =  (qltq2)  e  Q»  z  (•  i  q)  =  (w(- ;  q),  v(- ;  q))  is  the  solution  to  the  initial  value 
problem  (2.13),  (2.14)  or  (2.10)  -  (2.12)  with  El  set  equal  to  qt,  and  CqI  set  equal  to  q2. 

It  is  immediately  clear  that  the  optimization  problem  given  above  is  inherently  infinite 
dimensional.  The  admissible  parameter  set  Q  is  a  subset  of  a  function  space  and  die  evaluation  (and 
therefore  minimization)  of  the  least-squares  performance  index  |  requires  the  solution  of  an  infinite 
dimensional  evolution  equation.  The  introduction  of  finite  dimensional  approximations  is  essential 
to  the  development  of  practical  computational  methods.  Fundamental  to  the  approach  we  take  here 
is  a  weak,  distributional,  or  variational  formulation  of  the  initial  value  problem  (2. 13),  (2. 14).  We 
derive  the  weak  form  and  briefly  outline  existence,  uniqueness  and  regularity  results  for  solutions. 

In  the  usual  manner,  we  extend  the  operator  Cl(t)  given  by  (2.  IS)  to  an  operator  in  £  (V,  V) 
via 

(Cl(t)(v))(v) «  a  (t)(v,v),  v,v  e  V 
where  the  bilinear  form  a  (t)(-,-) :  V  x  V  -*  R  is  given  by 

a(O((ei.ij)>(02,X2))  =  <EI  L  Xi^O "  <EI  01*  lX2>0 “  <CDI  LXi-  l£2>o - 

(2.17) 

«*  f 

co(t,i)  J  e^xJdxfDj^i))  -  <o(t,)  J  e^xjdx,  x^  • 

o  o 

Standard  estimates  can  be  used  to  demonstrate  the  existence  of  positive  constants  k,  X  and  (3  for 
which 

kzftXv  j ,v2)l  £  k  Ivjly  lv2l«y  ,  vj  £  V,  i  *  1,2, 

and 

fl(t)(v,  v)  +  X|v|$t*|j|v|^/,  V  £  V,  t  E  [0,T1  . 

Consequently  (see  [27  ])  the  system  (2.13),  (2.14)  interpreted  as  an  initial  value  problem  in  V  or 
equivalently,  written  in  weak  form  as 


(2. 1 8)  <zt(t),v>H  *  a(t)(z(t),v)  +  <& (t),v>H 


v  e  V,  t  e  (O.TJ 


(2.19)  z(O)  =  zq 

admits  a  unique  solution  z  with  z(t)  e  V,  t  e  (0,T]  and 

z  e  L2((0,T);  V)  n  C([0,T1;  %)  n  H^O.T);  V). 

If  z  =  (w,v)  is  the  unique  solution  to  (2.18),  (2.19)  then 

u(t) «  L_1w(t),  t  e  [O.T] 

is  a  weak  solution  to  (2.6),  (2.7)  and  it  is  unique. 

Under  somewhat  stronger  hypotheses  than  those  given  in  Aj  and  A3  above,  the  existence  of 
strong  solutions  can  be  established.  Indeed,  if  in  addition  to  At  and  A4,  we  assume 

A2  The  mapping  t  — » a(t,  )  is  an  element  in  C*([0,T];  H^O.i))  for  some  T  >  0 

A3  The  mapping  t  -*  f(t,  )  is  an  element  in  C'tfO.Tl;  1^(0,  £))  and  g.h  e  C'fO.T]  (in  fact. 
Holder  continuity  will  suffice,  see  (291,  (37]) 

then  the  family  of  operators  {CUOl^o.T]  JPven  by  (213)  generates  a  unique  evolution  system 
(U(t,s) :  0  £  s  £  t  £  T)  on  3C  and  z  given  by 

t 

(2.20)  z(t)  «  Ud.Oz,,  +  J  U(t,s)3'(s)ds,  0  £  t  5  T, 

o 

is  the  unique  solution  to  the  inital  value  problem  (2. 1 3),  (2. 14)  and  satisfies  z(t)  e  49,  t  e  (O.T] 
with  z  e  C((0,T];  %)  n  C^((0,T];  X) .  Once  again,  with  z  *  (w,v)  now  given  by  (2.20), 
u(t) »  L_1w(t),  t  e  [O.T],  is  a  strong  solution  to  (2.6),  (2.7)  and  it  is  unique. 
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3.  An  Abstract  Approximation  Framework 

We  turn  next  to  •  discussion  of  a  general  approximation  framework  and  convergence  theory  for 
the  identification  problem  (ID)  formulated  above.  In  the  following  section  we  formulate  a  specific 
spline-based  scheme  to  which  the  general  theory  developed  here  applies. 

Fundamental  to  our  approach  is  the  construction  of  a  sequence  of  finite  dimensional  (with  regard 
to  both  die  state  dynamics  and  the  admissible  parameter  set)  approximating  identification  problems 
each  of  which,  under  appropriate  hypotheses,  can  be  shown  to  have  a  solution  that  in  some  sense 
(specifically,  subsequential  convergence)  approximates  a  solution  to  the  original  infinite 
dimensional  estimation  problem  (ID). 

In  the  discussion  to  follow,  we  exhibit  the  explicit  dependence  on  q  *  (qltq2)  €  <1  of  the  % 
inner  product  <  v>k  and  the  bilinear  form  d(t)(v)  given  in  (2.9)  and  (2. 17)  respectively  by 
using  the  notation  <-,->q  and  a(t;qXv).  For  each  Nj  -  1,2,...  and  each  N2  =  1,2,...  let 
Ni  N2 

W  and  V  be  finite  dimensional  subspaces  of  L^O.i)  and  V  respectively.  If,  for 

N  ">  “•  N 

N  *  (Nj,N2)  we  define  V  «  W  x  V  ,  then  V  is  a  finite  dimensional  subspace  of  both 
%  and  V.  Let  F  N  :  X  -*  V  N  denote  the  projection  map  of  X  onto  V  N  given  by 

(3.1)  V  (w,v)  *  (Pjw,  P^v  ) 

where  is  the  orthogonal  projection  of  1^(0, 2)  onto  W  1  and  is  the  orthogonal 

N. 

projection  of  H  onto  V  ,  both  computed  with  respect  to  the  standard  (unweighted)  inner 
products  on  the  respective  spaces  1^(0,  i)  and  H. 

The  Galerkin  equations  in  VN  corresponding  to  the  system  (2.18),  (2.19)  and  q  e<)  are 

(3.2)  <z*(t),  vN>  «  a(t ;  q)(zN(t),vN)  ♦  <& (t),  vN>  ,  vN  e  ,  0  £  t  £  T 


1  2 

For  each  Mj  e  Z+  =  { 1,2,...},  i  =  1,2,  let  SM  and  SM  be  finite  dimensional  subspaces  of 

*  S  2  1  2  1  2 
C[0,-t]  and  for  M  =  (Mj,M2)  define  by  <\M=  SM  x  SM  Let  0M  and  0M 

1  2 


denote  mappings  from  C[0,2]  onto  SMj  and  SMj  respectively  and  define  0^,  a  mapping  from 
3  onto  by 


1  2 

dM<q)  =  =  (dM^l)’  Q  =  (Ql*  Q2>  e  ^  • 


We  define  a  sequence  of  approximating  admissible  parameter  spaces  { QM} ,  Me  Z+  x  Z+  by 


Qm  =  dM*Q) 


and  formulate  the  sequence  of  approximating  identification  problems  as  follows: 


(ID^)  Given  £(t.)  e  Z4  =  1,2 . v,  find  (q^)  e  QM  which  minimizes 


JN(q)«  £  lr<zN(ti;q»-C(t)l’ 


over  Qm,  where  z^(-;  q)  is  the  solution  to  the  initial  value  problem  (3.2),  (3.3)  in  V^. 


N  K?  *N  *2  i  lJl* 

We  choose  bases  (8;  ,  { X;  and  for  the  finite  dimensional 


Nj  ^2  1  2  112  2 

spaces  W  ,  V  ,  SM  and  SM  respectively.  Then  qj^  e  SM  ,  qM  e  SM  and  ihe  solution 
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zN(  • ;  qM)  to  the  initial  value  problem  (3.2),  (3.3)  with  q  *  Qm  *  (Qm  •  Qm  )  can  be  written  as 


°m^m  • 


2  ^  , 


W;*  ,*^7,* 


and 


zN(t ; qM) * ; arf  +  X  ^  V! ; “m)^ •  t  e[0,T], 

i»l  i«l  i+Kj 


respectively.  Moreover,  ZN(- ;  aM)  is  the  solution  to  the  initial  value  problem  in 


«?♦«? 

R  given  by 


(3.5)  TIlN(aM)ZN(t)  *  AN(t ;  a*,)  ZN(t)  +  F^t),  t  e  (0,T] 

(3.6)  Z^OJ.zJ1. 


(fere  the  positive  definite  matrix  5tlN(aM)  is  of  die  form 


*V<V 


N 

7H  (o^) 


where  7n.^N(aM)  is  a  K* -square  matrix  with  components 


and  TTLjis  K^- square  matrix  with  entries 

ropy «  <m0£r.  xf>H 


For  each  t  £  0  the  matrix  AN(t ;  aM)  is  given  by  AN(t ;  aM)  =  AN(aM)  +  BN(t)  with 


y-  *  i 1  * ? 


E >M> 


An((*m)  = 


rEV  -c%m> 


BN(t)  = 


DN(t) 


where  EN(aM)  is  a  x  matrix  with  components 


[EN(oM)]i 


CN(aM)  is  a  -square  matrix  with  components 


[^(Om)] 


^  k+LM  k  _2  N  N 
ij  =  2Lr  “M  ^mDXj.DXj^ 


and  D^t)  is  a  x  matrix  with  components 


A 

[DN(t)]ij  =  cc(t,i)  J  e^xJdxDxffi)  -  <a(t ,  ■)  |  eJV) dx,  Dx*>0 . 


The  nonhomogeneous  term  FN(t)  is  given  by  F^(t)  =  (0,  F^(t))  where  F^(t)  is  a  vector  with 


entries 


[F?(t)li  =  <y0(o,  £><>„ . 


The  initial  data  is  of  the  form 


^-(oVc^.z^)1 


where  the  vector  Z^j  and  the  vector  zJJ,  are  given  component-wise  by 


<ii  =  <&%  ef>o 


-  <¥.  xf>H 


respectively,  and 


GN  = 


with  Gj*  a  K^-square  matrix  defined  by 

[Gfiij=<ef‘,ejl>0, 

and  G^  a  K^-square  matrix  with  components 

[G?]ij  =  <xfl.X?>H- 


It  is  now  easily  seen  that  the  finite  dimensional  identification  problem  (ID^,)  in  fact 
involves  simply  the  minimization  of  a  least-squares  performance  index  over  a  subset  of 

i1  »2 

+  N 

R  .  Furthermore,  the  evaluation  of  the  functional  J)  requires  only  the  solution 

N  N 

of  the  Kj  +  dimensional,  linear,  non-autonomous  ordinary  differential  equation  (3.5) 
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with  initial  conditions  specified  in  (3.6).  If  the  existence  of  solutions  to  the  finite  dimensional 
optimization  problems  can  be  established,  it  is  immediately  clear  that  they  can,  in  principle,  be 
computed  using  standard  techniques.  Conditions  which  guarantee  the  existence  of  solutions  to 
problem  (IDM)  and  the  fact  that  they  in  some  sense  approximate  solutions  to  the  original 
infinite  dimensional  estimation  problem  (ID)  are  given  in  the  following  theorem. 

Theorem  3.1  Suppose 

Hj  the  mappings  dM  are  continuous  from  Q  into  <\, 

Hj  for  each  q  e  Q,  dM(q)  — »  q  as  I M  I  — » «» with  the  convergence  being  uniform  in  q  for  q  e  Q, 
H3  the  spaces  VN  and  projections  are  such  that  if  {qN}  is  a  sequence  in  Q  with 

qN  -»  q  =  (qj,  e  Q  as  I N I  ->  «>  then  zN(t ;  qN)  -4  z(t ;  q)  in  1^(0, £)  x  H  for  each  t  e 
[0,T]  as  I N  |  — >  «>o  where  zN(- ;  qN)  is  the  solution  to  the  initial  value  problem  (3.2),  (3.3) 
with  q  =  qN  and  z(- ;  q)  is  die  solution  to  the  initial  value  problem  (2.18),  (2.19) 

corresponding  to  El  =  q}  and  CqI  =  q2. 

N  N  *  m  e 

Then,  each  of  the  problems  (IDM)  has  a  solution  (qM)  .  Furthermore,  the  sequence  {(qM)  } 

admits  a  ^-convergent  subsequence  whose  limit  q*  is  a  point  in  Q  and  is  a  solution  to  problem 
(ID). 

In  the  statement  of  the  theorem,  for  an  element  K  =  (Kj.Kj)  e  Z+xZ+  we  have  adopted  the 
notation  I K I  — » <»  to  denote  Kj,!^  — »  00 .  We  remark  that  it  is  also  true  that  the  limit  point  of 

any  ^-convergent  subsequence  {(q ^)  }  of  {(q^)*}  with  |  M  | ,  |  N  |  -» »  as  jjc  -> »  is  a 

solution  to  problem  (ID)  as  well.  Moreover,  if  problem  (ID)  has  a  unique  solution,  q*,  then  the 
N  *  * 

sequence  {(qM)  }  itself  converges  to  q  .  It  is  also  important  to  note  that  the  hypotheses  of  the 
theorem  do  not  require  that  Qm  cQ. 

We  have  established  results  analogous  to  those  given  in  Theorem  3.1  for  inverse  problems 
involving  parabolic  and  hyperbolic  systems  (see,  for  example  [12],  [13],  [16])  as  well  as  for 
related  methods  for  higher  order  equations  for  elastic  structures  (see  [4],  [5],  [6]).  For  the  flexible 
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structure  problems  treated  here,  the  essential  features  of  the  argument  remain,  for  the  most  part, 
unchanged.  We  therefore  only  briefly  sketch  them  below. 

Standard  continuous  dependence  results  for  linear  ordinary  differential  systems,  the  continuity 
assumptions  on  SJM  and  T  (and  therefore  on  dN  as  well)  and  the  fact  that  Q  is  a  compact  subset  of 
are  sufficent  to  conclude  that  there  exists  a  solution  (q^)*  e  QM  to  problem  (ED^) . 

The  definition  of  the  space  QM  (see  (3.4))  implies  the  existence  of  a  e  Q  for  which 

N  *  — N  — _m 

=  4M(qM)-  Since  Q  is  compact,  there  exists  a  subsequence  { q  .}  of  { qM }  with 

MJ  M 

— N*  *  -N^ 

q  j  — >  q  eQasjJc— >«o.  The  subsequence  { q  . }  can  always  be  chosen  with 
Kr  Kr 

ImJ  I,  I Nk  |  — »oo  asj,k— >oo.  It  follows  that 


n"  n“ 

D  «q^)*d  (q), 


qeQMi 


and  consequently  that 


*  N* 

U  ((q  -md  (d  .(q)). 


Assumption  H^,  above  and 


qeQ. 


l(q"V*q*l<v*  1^  j(q^)-q^l^+  Iq^-q’l, 

Nr  H  M1  S  MJ 


imply  (q  ;)  ->  q  as  j,k  ->  oo .  Taking  the  limit  as  j,k  ->  «>  in  (3.7)  above  with  an 

application  of  assumption  H3,  we  find  d(q*)  £  d(qX  q  e  Q ,  and  hence  that  q*  is  a  solution  to 
problem  (ID). 


In  this  section  we  outline  a  scheme  which  uses  piecewise  polynomial  spline  functions  and  show 


that  it  satisfies  the  conditions  and  hypotheses  of  Theorem  3.1.  We  first  treat  the  discretization  of 
the  admissible  parameter  set  Q. 

For  each  M  =  (Mj,  Mj)  eZ+xZ+  let  A*m  and  A^  denote  the  uniform  partitions  of  the  interval 
[0,-fc]  determined  by  the  meshes  [0,-t/Mj,  2-t/Mj, . . .  ,1 }  and  [O.-C/Mj,  21/ M^. . .  ,i } 
respectively.  For  m  =  1,2,. . .  and  A  a  partition  of  [0,  l]  let  Sp(m,A)  denote  the  usual  spline 

space  of  functions  in  C^m‘^[0,-t]  which  are  polynomials  of  degree  2m- 1  on  each  subinterval  of 

i  i  i 

A  (see  [36]).  We  then  define  SM  =  Sp(l,  A),  i  =  1,2.  In  this  case  we  have  dim  SM  =  Lj^  =  Mj  + 

i 

1,  i  =  1,2,  with  the  usual  "hat"  functions  forming  a  cardinal  basis  for  each  of  the  spaces  SM  ,  i  = 

i  i  ' 

1,2.  For  i  =  1,2,  let  S)M:  C[0,i]  — »  SM)  be  the  interpolation  operator  defined  by 

«.t)(4)-i(&)  •  j'0-1’2 . Mi 

t  1 

for  y  e  C[0,£].  The  theory  of  interpolatory  splines  (see  [33])  yields  the  continuous  dependence 
result 

I  "  Vfe  L  *  I  Yi "  Y2  L .  >  =  U 

where  Yj,  y2  e  C[0,£]  and  consequendy  that  hypothesis  Hj  of  Theorem  3.1  is  satisfied.  Also, 
the  approximation  result  (see  [36]) 

I  Vr-rloo  ^  w(y,  1/Mj) 

where  (o(y,8)  is  the  usual  modulus  of  continuity  of  ye  C[0,£]  with  respect  to  5,  together  with  the 
assumption  that  Q  is  a  compact  subset  of  =  C[0,-6]  x  C[0,-t]  and  the  Arzela-Ascoli  theorem 
yield  that  hypothesis  H2  is  satisfied  as  well. 

Next  we  define  a  state  approximation  and  verify  that  hypothesis  H3  holds.  As  above,  given 
N  =  (Nj,  N2)  e  Z+  x  Z+ ,  we  define  the  uniform  partitions  A^of  the  interval  [0,i]  determined  by 
the  meshes  [0,-6/Nj,  2-fc/Nj,. .  .,t },  i  =  1,2.  We  may  then  choose  either 

Nj  N 

W  1  =  Sp(l,  Aj ) 

N1  N 

W  =  Sp(2 4, )  • 


or 


In  the  first  case,  once  again  the  "hat”  functions  may  be  chosen  as  a  basis  with 


dim  W  =  Kj  =  Nj  +  1.  In  the  second,  the  standard  cubic  B-splines  (see  [33]), 


Nt  N,  + 1 

■>  M  * 


(Bj  )j  =  .j  ,  corresponding  to  the  partition  Aj  form  an  appropriate  finite  element  basis  with 


Ni 

dim  W  »  Nj  +  3.  In  either  case,  approximation  results  for  interpolatory  splines  can  be  used 


to  obtain 


(4.i)  Ippe-elo-^OasNj 


for  0e  1^(0,  t). 


We  set 


V  2  =  {(X(A  Dx(i),X)  e  H :  x  e  Sp(2,A  N ),  x(0)  =  Dx(0)  =  0} 


Then  V  c  V  and  defining 


N  N,  N,  N, 


Pi  =  B0  -  2Bj  -  2B  i  , 


*2  *2 


p.  =B.\i  =  2,3 . N2+l 


*n9  n,  n,  n. 

$.  =  (P.  2(i),Dp.  \X), p.  2),  i  =  1 ,2 . N,  +  1, 


,N2,N2  +  1 


the  collection  { (5.  }.=J  forms  a  basis  for  V  2  with  dim  V  2  =  K^  =  N2+1.  With  V 


as  defined  above,  it  is  not  difficult  to  show  (using  arguments  similar  to  those  in  [31]) 


Ip^&XJ-OI&XJIh-*0  as  N2 


,  V  '  A.'  «Xk>  .•>  .V  wN  .V  .N  A  V.V.\  .V.VA  .V 


for  Cn,£,x)  e  H  and 


(4.3)  IlP^X-LxIh-^0  as N2  -> oo 


for  £  e  V. 

So  as  to  avoid  obscuring  the  essential  features  of  our  argument  with  technical  details,  we  verify 
hypothesis  H3  for  the  spline-based  scheme  described  above  in  the  case  0  =  0.  The  term  which 
results  from  axial  loading  is  a  bounded  perturbation  and  does  not  involve  the  unknown  parameters. 
Showing  that  the  desired  convergence  continues  to  hold  in  the  presence  of  a  non-zero  axially 
directed  acceleration  requires  only  a  routine  extension  of  the  proof  which  we  give  below  (see  [14]). 

Suppose  that  {qN}  is  a  sequence  in  Q  with  qN  — >  q  e  Q  as  I N  I  — » °° .  Let  zN  =  (wN,  v*4) 
denote  the  solution  to  (3.2),  (3.3)  with  q  =  qN  and  let  z  =  (w,v)  denote  the  solution  to  (2.18), 
(2.19)  corresponding  to  q.  We  shall  require  the  assumption  that  z  is  a  strong  solution. 

In  the  estimates  which  follow,  we  simplify  our  notation  by  referring  to  the  inner  products 
(norms)  <•,•>  ^  ( M  N)  and  <v>-  ( M  >)onX  by  <v>N  ( M  N)  and  <•,•> 

(II)  respectively.  Also  note  that  with  a  =  0,  we  have  a  (t;  q)(v)  =  a  (q)(v). 

Since 

llzN  -  zll  £  llzN  -  ?PNzll  +  ll(PN  -  I)zll 

where  INI  denotes  the  usual  (unweighted)  product  norm  on  %  =  1^(0,  l )  x  H,  (3.1),  (4.1)  and 
(4.2)  imply  that  we  need  only  to  consider  the  term  llzN  -  PNzl! .  Letting  yN(t)  =  zN(t)  -  lPNz(t), 
using  (2.18),  (2.19),  (3.2),  (3.3)  and  the  fact  that  VN  c  V  we  find 

N  N  ,  _N  .  N  N  N 

<yt  .V  >N  -  <(z  -  V  z\  ,  V  >N  +  <Zt  ,  V  >  -  <zt ,  V  >N 

+  fl(qN)(yN,vN)  -  a  (qN)(z  -  1pNZivN)  +  a(qN)(z,vN)  -  a(q)(z,vN) 

+  <ST,vN>N  -  <CT,vN>  vN  e  VN,  0  <  t  £  T 

(4.4)  yN(0)  =  0. 


23 


Choosing  v**  =  y**  e  VN,  we  obtain 


=  <(z  -  FN*)t ,  yN>N  +  <(q!  -  q{Vt ,  w*  -  P}*w>0 
-  <q"L(v  -  Pjv),  wN  -  P^w>0  +  <q[*(w  -  pNw),  Lfv"  -  Pjv)>0  - 
<q^L(v  -  Pgv),  Lfv*4  -  P£v)>0  +  <(q»N  -  qi)Lv,  vF  -  P,Nw>0  - 
<(q"  -  qi)w,  L(v*  -  pjv)>0  -  <(q£  -  q2)Lv^(vN  -  Pjv)>0 . 

Recalling  that  Q  is  a  compact  subset  of  <\  and  that  for  q  =  (qlt  q2)  e  Q,  El  =  qt  and  CDI  =  q2  arc 
assumed  to  satisfy  assumption  Aj  of  Section  2,  we  find 

lllyVt  lu^-  ^IJs^lia-  A/ 

+  llyN||2+  lq,-<f  I*  Iwjj  +  |wN  -  P*w|j  + 

lLa-pJ)vlJ+lwN.pJ,w|j  +  ^la-p^wlj 
+  e  iLfv^  -  F^v)|j  +  -f-  |L(I- P^)v|j  + 

0  4e  4  u 

e|L(vN- P^v) |q+  IqJ^-qj  |2  |Lv|q  + 

|wN-pJw|j  +  .l  lq"-qj^  Mj  + 

eluCN-l£v)|J  +  ^|q"-qJ2|Lv|J 

+  e|L(vN-P^v)lj} 

where  Kq  is  a  positive  constant  and  e  is  an  arbitrary  positive  constant  Gathering  up  like  terms 
and  choosing  e  <-^  we  obtain 


1  IlyV  S  k,  {  11(1  -  FN)zlll2  ♦  I  ^  M  w  J  J  +  I  «I  -  pJ)  v  |  J  + 


I  (i  -  p^*>w  I  o +  itf-sJi  lL*lJ+  K,-51l!>l;+  UJ-  ^2'i|L;l2o>+K2NyN|12 

where  K,  and  iCj  are  positive  constants.  Integrating  both  sides  of  the  above  inequality  from  0  to  t 
and  recalling  (4.4)  we  obtain 

t 

(4.5)  ||yN(t)ll2  £  8  +  k2  J  ||yN(s)ll2ds 

o 

where 

T 

8  =  Kj  J* ( 11(1  -  F  )zg(s)ll2  +  I  q*  -  q,  I  *  I  wg(s)  I J  +  I UI  -  P2)v(s)  I J  +  1 0  -  P^)w(s)  I J  + 
o 

lq"-  q,  I*  |Lv(s)|j+  IqJ1-  q,  1^  |W(S)|2+  I  q2  -  q2 1  *  I  Lv(s)  I  J}ds  . 

Since  qN  -»  q  as  I N  I  — »  •»  and  z  ■ (w,v)T  was  assumed  to  be  a  strong  solution,  (4.1),  (4.2), 
(4.3)  and  (4.5)  together  with  an  application  of  the  Gronwall  inequality  yield  the  desired  result, 
llyN(t)ll— *0. 

A  close  inspection  of  the  estimates  above  reveals  that  they  depend,  to  a  large  extent,  on  the 
presence  of  the  viscous  damping  term  <CqI  l£j.  l£2>0  in  the  bilinear  form  a  (t)( , )  given  in 
(2. 17).  That  is,  we  require  that  q2  2:  a  >  0  for  some  a  >  0  for  all  q»(q,,  q2)  e  Q.  In  the 
absence  of  the  Voigt-Kelvin  damping  we  can  still  argue  the  convergence  of  zN  to  z;  however,  we 
must  assume  that  Q  is  H2-compacL  If  one  is  to  enforce  the  compactness  constraint  on  Q  when 
solving  the  finite  dimensional  optimization  problems  (a  desirable  implementation  feature  in  many 
cases  -  see  ( 10],[  1 1  ]),  this  stronger  assumption  becomes  especially  unappealing.  On  the  other 
hand,  by  employing  a  somewhat  different  (but  closely  related)  factorization  of  the  stiffness 
operator  Kq  than  the  erne  which  was  used  here  (one  which  is  formally  equivalent  to  rewriting  the 
initial-boundary  value  problem  (2.1)  -  (2.5)  as  a  first  order  system  in  the  states  El  D2u  and  ut  as 
opposed  to  and  ut)  hypothesis  H3  of  Theorem  3.1  can  be  verified  for  the  resulting  spline- based 
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scheme  under  the  present  assumptions  on  Q.  Unfortunately  this  scheme  is  also  difficult  to 
implement  and  from  a  numerical  standpoint,  has  not  performed  as  satisfactorily  as  the  one  based  on 
the  formulation  given  in  this  paper.  The  present  scheme  performed  well  whether  or  not  damping 
was  present  in  the  equation  and  hence  the  assumption  that  CDI  £  a  >  0  may  be  an  artifact  of  our 
proof  of  convergence  (see  Example  5.3  below). 

5.  Numerical  Findings 

We  present  and  discuss  some  of  the  results  which  we  obtained  from  our  numerical  studies  of  the 
scheme  that  was  described  in  Section  4.  All  codes  were  written  in  FORTRAN,  and  tested  and  run 
on  the  IBM  3081  at  either  Brown  University  or  the  University  of  Southern  California.  The  same 
codes  were,  with  only  minor  modification,  run  on  the  Cray  1  -S  at  Boeing  Computer  Services  in 
Seattle  with  support  made  available  to  us  through  the  National  Science  Foundation's  Super 
Computer  Initiative  program.  Examples  were  benchmarked  so  that  the  potential  benefits  of 
vectorization  to  our  research  program  could  be  accurately  and  effectively  assessed.  Our  findings 
are  described  below.  This  information  will  become  especially  important  to  us  when  we  begin  to 
consider  the  extension  of  our  general  approach  to  inverse  problems  involving  the  vibration  of  two 
dimensional  structures,  such  as  flexible  plates  or  platforms,  or  vibrations  of  structures  in  which 
nonlinearities  play  a  significant  role.  The  finite  dimensional  optimization 
problems  (IDM)  were  solved  using  the  IMSL  routine  ZXSSQ,  an  implementation  of  the  iterative 

Levenberg-Marquardt  quasi-Newton  algorithm.  The  finite  dimensional  initial  value  problems 
(3.5),  (3.6)  were  solved  in  each  iteration  of  the  minimization  procedure  (for  the  evaluation  of  the 
least-squares  performance  index  0  and  its  gradient  with  respect  to  the  parameters)  using  Gear's 
method  for  stiff  systems  (IMSL  routine  DGEAR). 

Our  codes  were  written  to  take  full  advantage  of  the  banded  structure  of  the  generalized  mass, 
stiffness  and  damping  matrices  afforded  by  the  use  of  polynomial  B-spline  elements.  All  necessary 
inner  products  were  computed  using  a  two  point  composite  Gauss- Legendre  quadrature  scheme. 

All  of  the  examples  presented  here  involve  fits  based  upon  displacement  measurements  obtained 
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through  simulation.  "True"  values  (which,  in  the  examples  below  will  be  denoted  with  an  asterisk, 
for  example  El*,  CpI*.  etc.)  for  the  unknown  parameters  were  chosen.  The  resulting  initial 
boundary  value  problem  (2.1)  -  (2.S)  was  then  solved  using  an  independent  integration  scheme. 
(We  used  a  seven  element,  quin  tic  spline  based  Galerkin  method  applied  directly  to  the  second 
order  system  (2.6),  (2.7)).  This  procedure  produced  sufficient  noise  in  the  data  so  that  the  use  of  a 
random  noise  generator  was  not  required. 

In  addition  to  the  test  example  numerical  studies  we  report  on  here  we  have  successfully  used 
methods  similar  to  those  developed  above  with  experimental  data.  These  results  are  presented  in 
detail  in  [9]. 

In  the  examples  which  follow  we  took  die  axial  loading  to  be  induced  by  an  acceleration  of  the 

base  or  root  of  the  structure  in  the  positive  x-direction.  In  this  case  we  have  (see  [34]) 

A 

0(ux) « -a^t)  [m  *  J  p(y)dy) 
o 

where  m  is  the  mass  of  the  tip  body,  p  is  the  linear  mass  density  of  the  beam  and  a,,  is  the  time 
dependent  base  acceleration. 

In  Examples  5.1  thru  5.4  below  we  took  l  «  1,  p(x)  *  3  -  x  for  0  £  x  £  1,  f(t,x)  *  exsin  2xt, 
g(t)  *  2e_t,  h(t)  *  e'^*,  aQ<t)  ®  1  for  0  £  t  £  1 .5,  8q0)  =  0  for  t  >  1 .5,  m  =  1 .5,  c  =  .  1  and  J  =  .52 
and  considered  the  estimation  of  the  flexural  stiffness  coefficient  El  and/or  the  viscoelastic  damping 
coefficient  CDI  only.  In  Example  5.1,  5.2  and  5.4,  the  fits  we  describe  are  based  upon 
observations  at  times  tj  =  .2i,  i  *  1,2, ...  ,5  at  locations  xj  *  .5,  .75  and  1.  In  Example  5.3 

observations  at  times  tj  =  .5i,  i  =  1,2,...,  10  at  locations  xj  =  .75  and  1  were  used.  In  all  of  the 

N,  N 

examples  we  discuss  here  the  space  W  was  generated  by  cubic  splines  (i.e.  as  Sp  (2.A, )) 

with  Nj  =  N2  =  N.  This  corresponds  to  the  approximation  of  the  first  and  second  components  of 

z  with  respectively  N  +  3  and  N  ♦  1  piecewise  cubic  C2  elements. 

The  compactness  constraints  on  the  spaces  Qm  were  not  enforced  when  the  finite  dimensional 

N 

optimization  problems  (IDM)  were  solved.  When  M,  and  became  large,  the  inherent 

ill-posedness  of  the  inverse  problem  became  apparent  as  the  performance  of  our  schemes 

deteriorated.  There  is  evidence  strongly  suggesting  that  this  situation  can  be  remedied  either  by 
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imposing  the  compactness  constraints  on  the  admissible  parameter  space  and  then  solving  the 
mnimiiMion  problem  using  a  constrained  optimization  procedure  (see  [  10),  ( 1 1])  or  by 
regularizing  the  least  squares  performance  index  (see  (24),  (25)).  We  intend  to  direct  our  attention 
to  these  ideas  in  the  near  future. 

EampltU 

In  this  example  we  consider  the  simultaneous  estimation  of  a  constant  flexural  stiffness 
coefficient,  El*  -  .15,  and  a  damping  coefficient  given  by  CDI*(x)  »y(1.5  -  tanh  (3x  -  1.5)), 
x  e  10,1),  with  y-  .01.  With  N  »  4,  M ,  ■  1  and  Mj  -  3  and  taking  start  up  values  (for  the  least 
squares  minimization  algorithm)  EI°  «  .  1  and  C^fx)  *  .015,  x  e  [0,1]  we  obtained  the  results 
shown  in  Figure  5. 1  below.  This  particular  run  required  approximately  30  seconds  of  CPU  time 
on  the  IBM  3081 


Figure  5. 1 

We  observed  that  how  well  the  scheme  performed  depended  upon  the  magnitude  of  the  scaling 
factor  y  As  y  was  decreased,  so  too  did  the  "sensitivity"  of  the  least  squres  performance  index  to 
the  damping  coefficient.  Results  similar  to  those  shown  in  the  figures  above  were  obtained  with 
Y  =  .005.  With  Y=  001,  on  the  other  hand,  we  were  unable  to  simultaneously  identify  both  of  the 
unknown  parameters.  However,  again  with  y  *  .001,  but  this  time  fixing  El  at  the  true  value,  we 


28 


were  able  to  identify  CqI  alone. 

*  •  1 

When  we  replaced  the  constant  El  with  the  linear  function  El  (x)»  1  -  jx  and  took 
Y*  1.  the  performance  of  the  scheme,  from  a  qualitative  point  of  view,  remained  unchanged. 


Example  5.2 

We  again  consider  die  simultaneous  estimation  of  the  stiffness  and  damping  coefficients.  We 
again  set  El*  *  .15  but  this  time  choose  CqI’Cx)  *  .01  (1.5  -  tanh  (20x  -  10)),  x  c  [0,1],  The 
identification  of  this  steeper  hyperbolic  tangent  function  has,  in  past  test  examples,  proven  to  be  a 
somewhat  stiffer  challange  for  our  methods  (see  [5),[6]).  With  N  *  4,  M,  *  1 ,  M2  =  3,  EI°  =  .  1 
and  CDl°(x)  *  .015  for  0  £  x  £  1,  we  obtained  the  estimates  which  are  plotted  along  with  the  true 
parameters  in  Figure  5.2. 


I  •>  •«  ••  •• 


Figure  5.2 

Also,  although  the  theory  was  not  explicitly  treated  here,  we  note  that  elements  other  than  linear 
splines  can  be  used  to  discretize  the  admissible  parameter  space.  Our  investigations  have  included 
numerical  studies  with  0-order  splines  (i.e.  piecewise  constant  functions)  and  cubic  spline 
elements.  Using  two  linear  elements  to  approximate  El  (i.e.  M,  =  1 )  and  nine  cubic  elements  to 
discretize  CD1  we  obtained  the  estimates  shown  in  Figure  5.3.  We  have  obtained  an  acceptable 


estimate  for  CpI  with  as  few  as  six  cubic  elements. 

In  the  tests  reported  on  for  the  present  example,  residuals  were  typically  in  the  range  10"^  to 
10-8  with  CPU  times  from  25  to  40  seconds. 


Figure  5.3 


Example  3J 

In  this  example  we  identify  only  the  spatially  varying  flexural  stiffness  coefficient 
EI*(x)  *  1.5  -  tanh  (3x  -  1.5),  x  e  [0,1],  in  a  model  with  no  viscoelastic  damping  (CDI  =  0).  In 
Section  4  we  remarked  that  our  convergence  arguments  required  either  the  presence  of  viscoelastic 
damping  in  the  model  or  that  the  admissible  parameter  set  Q  be  compact  in  the  stronger  H2 
topology.  The  results  shown  in  the  figure  below  suggest  that  this  is  only  an  artifact  of  our  proof 
and  not  a  fundamental  requirement  for  the  convergence  of  our  approximation  (i.e.  the  absence  of 
damping  does  not  appear  to  effect  the  overall  performance  of  our  scheme). 

Taking  N  •  4  and  M,  equal  to  1  thru  8  we  produced  the  results  shown  in  the  series  of  graphs  in 
Figure  5.4.  The  initial  estimate  or  start-up  value  for  El*  was  taken  to  be  the  constant  function 
EI°(x)  «  1  forOS  x£  1. 


Recalling  our  earlier  remarks,  the  oscillations  which  appear  in  the  graphs  corresponding  to 


Mj  =  6, 7  and  8  due  to  the  inherent  ill  posedness  of  the  estimation  problem  are  not  unexpected. 

In  fact,  as  or  M2  —»<*>,  the  appearance  of  the  undesirable  oscillations  in  our  final  estimates 
occurred  in  virtually  every  test  we  ran.  As  we  have  noted  earlier  however,  preliminary  findings  in 
related  studies  [10]  and  [1 1]  regarding  the  enforcing  of  the  compactness  constraints  and  the 
subsequent  use  of  constrained  optimization  techniques  to  solve  the  approximating  finite 
dimensional  identification  problems  suggest  that  this  difficulty  can  be  overcome.  Our 
investigations  in  these  directions  are  continuing. 

In  addition,  the  series  of  tests  corresponding  to  the  graphs  in  Figure  5.4  were  benchmarked  on 
the  IBM  3081  and  the  Cray  1-S.  The  same  estimates  were  obtained  on  both  machines.  However, 
we  were  able  to  achieve  a  speed-up  factor /,  of  7  - 10  on  the  vector  machine.  The  CPU  times  are 
reported  in  Table  5.1.  In  comparing  the  CPU  times  on  the  3081  for  this  example  with  the  times 
reported  for  the  previous  examples  it  is  important  to  note  that  the  results  here  were  based  upon 
observations  taken  over  the  longer  time  interval,  [0,5],  versus  the  interval  [0,1]  for  examples  5.1 
and  5.2. 

Example  5.4 

In  Figure  5.5  below  we  plot  the  final  estimates  obtained  when  we  attempted  to  use  our  scheme 
to  simultaneously  identify  the  spatially  varying  flexural  stiffness  coefficient 
EI*(x)  =  .5  +  4x(l  -  x),  x  e  [0,1],  and  viscoelestic  damping  coefficient, 

CDI*(x)  =  .1  (1.5  -  tanh  (3x  -  1.5)),  x  e  [0,1] .  The  start-up  values  for  the  iterative  least  squares 
minimization  routine  were  taken  to  be  the  constant  functions  EI°(x)  =  1  and  CDI°(x)  =  .15  for  0  £  x 
<;  1 .  The  graphs  in  the  figure  were  obtained  with  N  =  4  and  a  linear  spline  discretization  of  the 
admissible  parameter  set  Q  with  Mj  *  M2  =  3.  In  all  of  our  tests  with  this  example  the  minimum 
sum  of  the  squares  of  the  residuals  was  in  the  range  10'7  -  10‘8  with  the  optimization  typically 
requiring  50  -  70  seconds  of  CPU  time. 
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IBM  3081 
(CPU  sec.) 

CRAY  1-S 
(CPU  sec.) 

110 

12.5 

164.1 

20.8 

207 

23 

249 

32 

245.5 

36 

404.4 

41.6 

346.5 

45.6 

437.9 

49.1 

For  this  example  we  also  tried  a  cubic-spline  based  discretization  for  Q.  We  considered  all 

)|l  ]|c 

possible  combinations,  linear  splines  for  El  together  with  cubic  splines  for  CDI  ,  cubic  splines  for 

El*  together  with  linear  splines  for  CDI*,  etc.  Although  small  values  for  the  sum  of  the  squares  of 
the  residuals  were  obtained  in  each  instance,  our  by  far  best  approximation  to  the  true  parameters  is 
the  one  shown  in  Figure  5.5  which  corresponds  to  a  linear  spline  based  discretization  for  both 
components  of  the  admissible  parameter  set 

Holding  CDI  fixed  at  the  true  value  and  using  cubic  splines  to  identify  El  and  then  holding  El 

fixed  at  the  true  value  and  using  cubic  splines  to  identify  CDI*  we  were  able  to  obtain  the  estimates 
plotted  in  Figures  5.6  and  5.7  respectively.  The  estimate  for  El  graphed  in  Figure  5.6  was 

j|r 

obtained  with  10  cubic  elements  while  the  estimate  for  CDI  in  Figure  5.7  is  a  linear  combination 

of  6  cubic  elements.  An  inspection  of  these  figures  reveals  that  while  the  approximations  obtained 
are  at  least  marginally  acceptable,  it  is  also  not  surprising  that  our  scheme  had  some  difficulty  when 
we  attemped  to  identify  both  parameters  simultaneously  with  a  cubic  spline-based  discretization  for 
either  one  or  both  components  of  Q. 

For  this  example  we  also  looked  at  the  robustness  of  our  iterative  scheme  with  respect  to  the 

initial  values  chosen  (i.e.,  EI°  and  CpI0).  In  Figure  5.8  we  plot  those  points  in  the  Cpl°  -  EI° 
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plane  which  correspond  to  the  start  up  values  we  tried.  The  point  marked  with  "  *  "  corresponds  to 
the  startup  values  which  produced  the  approximations  shown  in  Figure  5.5.  The  points  marked 


with  "  x  "  correspond  to  start-up  values  which  led  to  essentially  the  same  extimates  as  those  shown 

in  the  figures.  The  points  marked  with  an "  0  "  correspond  to  start-up  values  for  which  the 
scheme  did  not  converge.  The  region  whose  boundary  is  denoted  with  dashed  lines  corresponds  to 
a  "convergence  envelope"  for  the  vector  valued  function  ( CpI  ,  El  ).  An  analogous  study  was 
carried  out  for  Example  5.2,  for  which  similar  robustness  results  were  obtained. 


Figure  5.8 


Finally  we  offer  several  summary  comments  on  some  of  our  other  numerical  findings.  In 
virtually  all  examples  we  tried,  we  found  that  the  estimates  yielded  by  the  scheme  which  we 

develop  here  based  on  state  space  coordinates  (D^j.Uj)  and  the  ones  yielded  by  the  scheme  based 

on  a  state  space  formulation  in  coordinates  (u,ut)  described  in  [5]  and  [6]  were  comparable. 

Although  in  any  given  example  one  scheme  or  the  other  may  produce  a  somewhat  better 

approximation  to  the  true  parameters,  we  found  it  impossible  to  designate  or  identify  a  clear  favorite 
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among  the  two  methods. 


We  also  ran  a  series  of  tests  in  which  we  varied  the  boundary  conditions  at  the  free  end  of  the 
beam.  That  is,  in  addition  to  the  tip  body  end  condition  we  considered  a  beam  which  is  clamped  at 
one  end  and  free  at  the  other  with  either  a  point  mass  (c  =  J  =  0)  or  no  mass  (m  =  c  =  J  =  0)  rigidly 
attached  at  the  tip.  We  also  studied  the  effect  that  the  presence  or  absence  of  external  forces  and/or 
moments  at  the  tip  of  the  beam  (i.e.  g  and  h)  has  on  the  performance  of  our  scheme.  Based  upon 
these  tests,  we  found  it  difficult  to  make  definitive  statements  regarding  "best"  experimental 
procedures  for  identification  of  structural  parameters  with  our  schemes.  However,  we  are  able  to 
offer  several  observations.  For  example,  with  a  point  mass  at  the  tip,  the  schemes  performance 

was  enhanced  when  an  external  moment  was  applied  at  the  tip  (i.e.  h  *  0).  On  the  other  hand,  the 

presence  of  an  externally  applied  force  in  the  transverse  direction  (i.e.  g  *  0)  did  not  appear  to  have 
any  effect  at  all.  Also,  with  no  mass  at  the  tip,  the  scheme  was  most  effective  when 
g  =  h  =  0.  In  general  we  found  the  scheme  to  be  most  dependable  with  tip  body  end  conditions. 
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